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I Abstract. The solar chromosphere is very dynamic, due to the presence of large amplitude hydro- 

dynamic waves. Their propagation is affected by NLTE radiative transport in strong spectral lines, 
which can in turn be used to diagnose the dynamics of the chromosphere. We give a basic intro- 
duction into the equations of NLTE radiation hydrodynamics and describe how they are solved in 
current numerical simulations. The comparison with observation shows that one-dimensional codes 
can describe strong brightenings quite well, but the overall chromospheric dynamics appears to be 
governed by three-dimensional shock propagation. 
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^ INTRODUCTION 

At low spatial and temporal resolution, the solar atmosphere has a well-defined average 
^ , temperature structure (Vernazza et al. [1], Fontenla et al. [2]; for snapshots of the 

\^ • most recent variants of such models see also [3, 4]): from the surface (i.e., the layer 

^ . seen in continuum radiation in the visible part of the spectrum) the temperature drops 



O 



continuously within the photosphere, from around 5 800 K down to a minimum value 
of the order of 4 000 K at a height of about 500 km. Beyond the minimum, the average 
^ ■ temperature increases first gently in the chromosphere, then rapidly in the transition 

l>- ! region, which starts at a height of some 2000 km, reaching nearly a million K even 

in the coolest parts of the corona, and several million K in the hottest parts. (For an 
extended review of the solar atmosphere see e.g. Solanki and Hammer [5].) This outward 
^ '. temperature rise must be supported by mechanical heating due to a combination of 

^ [ magnetic and nonmagnetic mechanisms (as reviewed e.g. by Narain and Ulmschneider 

[6, 7]; Ulmschneider and Musielak [8]), which are ultimately powered by the convective 
motions below the photosphere. 

With the advent of more advanced instrumentation, the achievable spatial and tem- 
poral resolution has been enhanced, and we have increasingly realized that this crude 
picture can serve only as a rough guide to the average solar atmosphere. The real atmo- 
sphere is characterized by a high level of fine structure, which in the upper chromosphere 
(and even more in the corona) is mostly due to thermal effects within magnetically con- 
fined plasma ([9, 10]). Moreover, the chromosphere is extremely dynamic. For example, 
much of the emission from the upper chromosphere comes from fine structures like 
spicules and fibrils (e.g., [11, 12, 13, 14]), in which the plasma is accelerated to high 
velocities, up to several times sound speed. 

Even outside of active regions, the solar surface is permeated by magnetic fields. 



Chromospheric Dynamics and Line Formation 



May 25, 2008 



1 



Stronger, long living magnetic flux concentrations tend to be arranged by the super- 
granular flow in a network-like pattern, which is visible up into the transition region. 
The interior parts of cells in this chromospheric network are not entirely field-free, but 
the magnetic field is weaker than in the network and very likely unimportant for the 
dynamics. If, and to what extent, it affects the heating of the plasma is still under dis- 
cussion. The importance of the magnetic field changes with height, since magnetic flux 
concentrations expand with height and ultimately fill all available space in the upper 
chromosphere and corona. The current paper will concentrate on the dynamics of these 
cell interior parts of the solar chromosphere; therefore we neglect the effects of the mag- 
netic field in the equations to be derived in subsequent sections. 

The connection between the solar chromosphere and the underlying photosphere has 
recently been studied by Woger [15], who also produced a movie clip [16] that shows the 
photo spheric and chromospheric dynamics both in the network and in the cell interior at 
high spatial resolution. 

Even those parts of the solar chromosphere where magnetic structuring is unimportant 
are known to be permeated by strong waves. They develop out of waves of relatively 
small amplitude that arise quite naturally in the underlying convective layer. When these 
waves travel upward, they experience a decrease in density p, due to the gravitational 
stratification. This leads to an increase of the wave amplitude Vmax, since the wave 
energy flux F^ave °= PVmax<^f where Cg is the sound speed, is conserved as long as the 
waves do not dissipate and if they are not damped by radiative energy exchange between 
wave peaks and valleys. The larger the wave amplitude becomes, the more important 
are nonlinear effects. As a result of the latter, the wave peaks propagate faster than the 
valleys and try to overtake them; thus the waves "break" and form quasi-discontinuous 
shocks, which dissipate the wave energy by thermal conduction and viscosity in their 
steep shock fronts (e.g., [17]). 

Most of the solar radiation in the visible originates from the photosphere; the chromo- 
sphere is virtually transparent at these wavelengths. For ground-based observations of 
the chromosphere one is restricted to the central regions of a few very strong absorption 
lines in the visible, near UV, and near IR; while from space one can observe the chro- 
mosphere also in EUV emission lines and continua. The transition region and corona, 
finally, emit predominantly in the EUV, X-ray and radio ranges. 

The most important spectral lines for studying chromospheric dynamics are the H- 
and K- lines as well as the infrared triplet lines of singly ionized calcium (Ca ll). The h- 
and k-lines of Mg ii would also be ideal for this purpose, but they lie too far in the UV to 
be observable from the ground. These strong lines are not only of high diagnostic value, 
but (along with numerous iron [mostly Fe II] and some hydrogen lines) they represent 
also the dominant cooling agents of the solar chromosphere (Anderson and Athay [18]) 
and must therefore be treated adequately in numerical simulations. 

Fig. 1 was derived from a time series of the Ca ii H line profile obtained in August 
2005 with the Echelle spectrograph of the German VTT telescope at the Observatorio 
de Tenerife [19]. As a result of the dynamic behavior of the chromosphere during 
the time series, the line profile is highly variable: Perturbations are seen all the time 
- sometimes small ones moving towards the core, but more commonly one notices 
large parts of a line wing change synchronously. A movie of the line profile during 
the observation demonstrates these changes (Rammacher et al. [20]). Fig. 1 shows the 
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FIGURE 1. Inner blue wing of the very strong Ca II H line. The middle curve (dotted) shows the time- 
averaged profile during a 4 300 s observation sequence. Each part of the Une profile probes different parts 
of the solar atmosphere, from the mid-chromosphere in the line core at 3968.49 A down to the lower 
photosphere in the far line wing (not shown). The indicated spectral range contains also a number of weak 
Fe I lines, which are formed in the mid-photosphere. When waves propagate through the atmosphere, they 
affect the intensity at the respective wavelengths. The upper and lower curve show the maximum and 
minimum values that were reached during the observation sequence. 



time- averaged profile and the minimum and maximum intensity values that occurred 
during the sequence. 

In these lectures, we introduce the physical principles that govern the highly dynamic 
behavior of solar and stellar chromospheres, thereby restricting ourselves to those re- 
gions where magnetic fields do not dominate the dynamics. These physical principles 
are described by the time-dependent equations of radiation hydrodynamics. To outline 
their theoretical foundation, we first derive the basic hydrodynamic and thermodynamic 
equations. Next we discuss elements of radiation theory, first the basic concepts and then 
the complications that arise in chromospheres. Finally we give a brief overview of how 
these equations are solved in numerical calculations and to what extent current state-of- 
the-art simulations can describe the observed chromospheric dynamics. All main chap- 
ters are preceded by recommendations of literature suitable for further reading. 



Chromospheric Dynamics and Line Formation 



May 25, 2008 



3 



BASIC HYDRODYNAMICS AND THERMODYNAMICS 

General references on the hydrodynamical, thermodynamical, and mathematical con- 
cepts treated in this chapter include [17, 21, 22, 23, 24]. 



Continuity Equation, Euler Frame 




FIGURE 2. Euler frame Eg (left): observe the mass flow out of a fixed volume V; Lagrange frame El 
(right): follow a mass element M 

In the Euler frame Le we consider a Cartesian coordinate system with origin O and 
unit vectors Cx, e^, e, in x-, y-, z-directions. We observe the density p(r,?) [gcm^^ ]^ 
gas pressure p{rj) [dyncm^^], and temperature T(rj) [K] of plasma flowing out of 
a fixed volume V with velocity v(r,?) [cms^^ ] (Fig. 2). r is the radius vector with 
components and t is the time. Let dA be a directed surface element of V. Then 

the amount of matter flowing out of V across its surface per unit time (Fig. 2) is 



^ pv ■ dA = Jv-p\dV 



(1) 



where the equal sign is due to Gauss's theorem. The amount of matter missing per unit 
time from volume V is given by 



dt 



pdV 



dt 



dV 



(2) 



where the equal sign is due to the fact that V does not depend on time. Equating Eqs. (1), 
(2), and because V is arbitrary one finds the continuity equation: 



dt 



(3) 



We consider it instructive to specify, in square brackets, units for newly introduced quantities. This 
proves to be particularly helpful for radiative terms introduced below. We choose CGS units. 
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Lagrange Frame 



In a Lagrange frame El we monitor the mass element M (Fig. 2) initially contained 
in V (at time t = and position r = a), which at a later time t is in volume V' at radius 
vector r(a,?) = (;c(a,?), j(a,?),z(a,?)). Vector a uniquely identifies amass element at the 
initial time t — 0. We consider physical variables such as p on an (x,t)-plane, on which 
the position of the mass element M describes a path x{t), where x{t — 0) = a, so that 
p = p{x{t),t). The rate of change of p along that path can be pictured as the sum of the 
incremental variations of p in and jc-directions, 

where the sign (=) means equal by definition and 

^' = ^ (5) 
at 

is the velocity of the mass element in ^-direction. The time derivative d/d? = {d /dt)^^ 
in the Lagrange frame is the so called substantial derivative, sometimes also called 
material or total derivative. It describes the time-dependence of a physical quantity in a 
moving mass element. In three dimensions, and for any physical function h (such as p, 
p, T, etc.), the substantial derivative is analogously given by 

dh 

- = ^ + v.V/, . (6) 



Equation of Motion 

In the Lagrange frame Z^, according to Newton's third law the inertial force (= mass 
dM times acceleration dv/d?) is balanced by the forces acting on the mass element, 

M V V V 

Here the acting forces consist of a volume force density f [ dyn cm^^ ] and a pressure 
force. In stellar atmospheres the volume force density is due to gravity, 

,= -^e. = -pje, , (8) 

where is the unit vector in radial direction, the stellar mass, G = 6.614 x 10^^ 
cm^ g^ s^ the gravitational constant, and r the radial distance from the center of the 
star. As the extent Ar of stellar atmospheres is often small (Ar <^ r), we can replace 
GM^/r^ by a constant gravitational acceleration g, e.g. g = 2.74 x 10"^ cms~^ for the 
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solar atmosphere. Using Eqs. (6), (7), and (8), and noting that the volume of the mass 
element is arbitrary, we obtain the equation of motion 



P (^^ + v-Vvj =-V/7-pge, . (9) 

In this equation we have neglected viscosity and radiation pressure. It can be shown that 
both are unimportant in the atmospheres of typical late-type stars (i.e., stars of spectral 
type beyond late A). 



Energy Equation 

It is convenient to picture the moving mass element as enclosed by a wall and 
to carefully monitor the energy flowing through that wall. A powerful book-keeping 
quantity for this is the specific (i.e., per unit mass) entropy S [ergg~^ K^^ ]. If the walls 
are impermeable, energy conservation just means S = constant in El. This can be written 
with Eq. (6) 

d5 dS ^ 

— = — +v-V5 = . (10) 
at dt 

In atmospheric regions close to the star (in photospheres and most of the chromospheres) 
only two processes by which mass elements gain energy are important, radiation and 
Joule heating. In coronae and some chromospheric regions one has steep temperature 
gradients and sometimes large velocity gradients. Here thermal conduction and viscosity 
are other important heating (or energy loss) mechanisms. One thus has the relation 



dS dS „^ dS 

:r = ^+v-V5= — 

dt dt dt 



= ^ + ^ + ^^ + ^X , (U) 
pT pT pT pT ' 



called entropy conservation equation. Here dS/dt\^^^ [ergg^^ K^^ s^^ ] is the heating 
function resulting from external heating by radiative, Joule, thermal conductive, and 
viscous heating, where Or, 4>c, and Oy [ergcm^ s^ ] are the net radiative, thermal 
conductive, and viscous heating rates, respectively. These functions will be discussed 
below. Joule heating, ^>j, will not be considered here, although magnetic fields and 
associated heating (magnetohydrodynamics, MHD) are important at least in some parts 
of stellar chromospheres, as discussed above. 



Atmospheric Gas Composition 

Stellar atmospheres consist of ideal gases, which are described by the universal gas 
law 

p = Y^nikT = p— , (12) 

i 

where k = 1.3807 x 10~^^ ergK~^ is the Boltzmann constant, /i [gmol~^ ] the mean 
molecular weight, and % = 8.3145 x 10^ ergK~^mol~^ the universal gas constant. 
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Hi are the number densities [cm^-^] of the different types of particles i (atoms, ions, 
and electrons). With a given mixture of chemical elements in the stellar gas the mean 
molecular weight can be computed as a function ix — ix{T,p). For solar type neutral gas 
jU 1.24, while for fully ionized gas jU ^ 0.60. 



Basic Elements of Thermodynamics 

In many applications in stellar photospheres, chromospheres, and coronae it is suffi- 
cient to consider atmospheric gases that are either neutral or fully ionized. Under these 
conditions the thermodynamic relations are particularly simple. The specific internal 
energy Eg [ erg g~^ ] is given by 

Es = c,T , (13) 
where Cy [ erg g~^ ] is the specific heat per unit mass for constant volume, given by 

1 9t 



r-iju ' 



(14) 



and y=Cpl Cy is the ratio of specific heats, which is constant (7 = 5/3) for either neutral 
or fully ionized gases. With the specific volume Vj [ cm-^ g~^ ], 

V,= ^ , (15) 
P 

and the fundamental laws of thermodynamics we obtain 

rd5 = d£', + /7dV, = d£',-^ dp . (16) 

From these equations the relationships between the thermodynamical variables 5, Eg, 
T, p, and p can be computed. Specification of two of these allows to determine the 
remaining variables. 



Conservation Equations 

It is often convenient to write the hydrodynamic equations in conservation form, in 
terms of conserved quantities / (mass, momentum, energy) and associated fluxes F 
(mass flux, momentum flux, energy flux): 

§^ = -V-F + C , (17) 

where C is a source term. 
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a. Conservation of mass 

The continuity equation (3) is already in conservation form 

dp 



V.pv. (18) 



b. Conservation of momentiim 

From Eqs. (3) and (9) one obtains 



dpy dy dp „ „ „ 

-^ = p^ + v^ = -pv-Vv-vV-pv-Vp-pge, , (19) 
which can be written 

^ = -V-(pvv + /7^)-pge, . (20) 
Here vv is a dyad (e.g., [23, 24]), and 

U = exex + eyey + ez^z , (21) 

is a unit dyad. 

c. Conservation of energy 

There are three types of energy densities [ erg cm~^ ] for gas elements in nonmagnetic 
stellar atmospheres: 

pEs internal energy (energy of microscopic undirected motion of atoms and ions), 
^pv^ kinetic energy (directed motion of the entire gas element), 
p<j> potential energy (gas element in the gravitational field). 
The gravitational potential is given by 

= *- . (22) 

r 

With this and Eq. (8) the volume force density can be written 

f=-pge, = -p^e, = -pV0 . (23) 

The time derivative of the total energy density can be written as a sum of four terms 

d f I r, \ 1 9(9p dv dpEs dpd) 

5j(2Pv^ + p£. + P*j=2V^^+pv.^ + ^ + ^ . (24) 

Modifying the terms on the RHS using Eqs. (3), (9), (11), (16), and (23), one obtains the 

energy conservation equation: 

I (ipv^ + p£. + p^) = -V.pv (iv^ + £, + ^ + ^) +pr ^ 
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= -V-pv(^-v2+£, + ^ + (^j +4>R + 4>j + ^>c + ^>v • (25) 

It is seen that there are three energy flux components [ erg cm~^ s~^ ]: 
pvjV^ kinetic energy flux, 
p V {Eg + p/p) enthalpy flux, 
pv0 potential energy flux. 



Heating by Viscosity and Thermal Conductivity 

The particle transport processes which occur in the presence of velocity and tempera- 
ture gradients contribute to the local heating. 
The thermal conductive heating rate 4>c [ erg cm~^ s~M is given by 

^c = -^Kjh^ . (26) 
ax ax 

The viscous heating rate ^>v [ erg cm~^ s~M is 

^Y = r],iJ^] . (27) 



\dx 

The coefficients of thermal conductivity f^/, and viscosity 77v,i are functions of T and p 



ELEMENTARY RADIATION THEORY 

Using the three time-dependent hydrodynamic equations (18), (20), and (25) together 
with the ideal gas law Eq. (12), the thermodynamic relations and functions such as 
IJ.{T,p), Kth{T,p) and r\yis{T,p), one would be able to compute the dynamics of stellar 
chromospheres, except that radiation is critically important. For this reason we now give 
a short review of the basic equations of radiation theory and derive the radiative heating 
rate. For further reading we suggest in particular [25, 26], but also [27, 28], and for 
looking up equations and atomic data [29, 30]. 



Basic Radiation Quantities 



a. Solid angle 

Consider an orthogonal Cartesian coordinate system [x, y, z). Let the latitude angle i? and 
azimuth angle (p be defined as seen in Fig. 3, left. The solid angle is a surface element 
on a sphere with unit radius 

da = sint?dt?d(p , (28) 
and has the dimension [ sr ] (for ste radian). 
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FIGURE 3. Left: SoUd angle. Right: Intensity 



b. Intensity 

The intensity is the basic quantity of radiation theory. Consider the energy of all 
photons that flow through a window of area AA in normal direction e„ into the solid 
angle AO. per time interval At and frequency interval Av (Fig. 3, right). Assume that r is 
the radius vector that points from the origin O of a coordinate system (with unit vectors 
Gx, Cj, Cz) to the location of the area AA. Then the monochromatic intensity ly, also called 
specific intensity, is given by 

The limit is taken for AA,A?,A^2,Av — > 0. /y has the dimension [ergcm^^s^^ sr^^ 
Hz-i ] = 10^3 [ Wm-2 sr-i Rz^ ]. RecaU that 1 J = 1 Ws = 1 Nm = 1 kgm^ s^^ = 
10^ erg. 

c. Mean intensity 

The mean intensity is the specific intensity averaged over all angles, 

Tt 2% 



Jv (r, t) = (^ly (r, e„,t) dQ. = -^J Jly sin-& dt? d^ . 



(30) 





In most cases there is no dependence of ly on (p, and with the angle cosine 

jl = 008-0- (31) 

the mean intensity can be written 



Jv{r,t) = ^Jlvsmjd^=^Jlyd^ . 



(32) 



-dp -1 

The dimension of Jy, [ erg cm^^ s^^ sr^^ Hz^^ ], is the same as that of ly. 
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FIGURE 4. Radiative flux 



d. Radiative flux 

Consider a window with area AA and normal e„ (Fig. 4). Assume that from both sides 
photons flow in arbitrary directions through AA. The net energy IS.E transported in e„ 
direction by all the photons flowing through AA per time interval A? and frequency band 
Av is the radiative flux given by 

^^^"''"''^ =^^Ja^ ' ^^^^ 

where the limit is taken for AA, A?, Av 0. Fy has the dimension [ erg cm^^ s^^ Hz^^ ] 
and can be derived from ly by summing over all projected contributions from the 
individual light rays in e', direction, 

K In 

Fy{r,en,t) = ^Iv{r,e'„,t)e'„-endQ' = J J ly cos sin ■& d-&^ dq) , (34) 

M -d^i 

where '& is measured from the e„ direction. Here one takes into account that photons 
from see only a projected area AAe^ ■ e„. In most cases there is no (p dependence, and 

+1 

Fy (r,e„,0 = 27t J ly II dji . (35) 
-1 

Note that Fy and ly are scalar quantities that depend on the normal vector e„ of the 
considered directed unit area. In an isotropic radiation field one has Fy = 0. 



Radiation Field, Level Populations, LTE and NLTE 

Consider the interior of a cavity which has been submerged for a long time in a gas 
of constant temperature T. All time-dependent processes have quieted down, and a final 
state has been reached, called thermal equilibrium (TE). It is described by four important 
relations: the Planck function, the Boltzmann distribution, the Saha equation, and the 
Maxwell velocity distribution. 
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FIGURE 5. Energy levels 

a. Planck function 

In TE the intensity is given by the Planck function: 



_ 2hv^ 1 

= ^ ^hv/kT _ 1 ' (36) 



which upon frequency integration gives 

Iydv=B= [Bydv = -T^ . (37) 

B is called the integrated Planck function. Here h = 6.626 x 10~^^ ergs is the Planck 
constant and o = 5.671 x 10^^ erg cm^^ s^^ K^^ the Stefan-Boltzmann constant. 

b. Boltzmann and Saha equations 

Consider the energy levels of the atoms or ions of a gas (Fig. 5). The number ni of 
atoms or ions per cm-' in the bound energy level / is called the population of level I. 
nj^ is the population of the continuum, i.e., the number of atoms or ions per cm^ having 
their electron removed by ionization. Bound levels in TE are described by the Boltzmann 
distribution: 

m gi gi 

where gi are statistical weights (for hydrogen e.g. gi = 2i^) and is the energy 
difference between the levels (for tabulated values of these quantities see Allen [29, 30]). 
Continuum levels in TE follow the Saha equation: 

3 /2 

nk ne ^2ukf 2Ttm^kT \ ^-Ei/kT ^^g^ 



m gi V 

where (see [29, 30]) Uk is the partition function, E] the ionization energy from level /, 
nip = 9.1094 X 10^28 gjj^g 

mass of an electron, and n^ the number of electrons per cm-^^. 
Note that the number of particles per cm^ is usually called the number density of these 
particles. 

c. Maxwell distribution 

In TE, atoms, ions, and electrons also obey the Maxwell velocity distribution: 

n \27tkT J 
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Here dn(v) is the number of particles per cm^ with mass m and velocities in the interval 
V to V + dv, while n is total number of these particles per cm^, irrespective of velocity. 

TABLE 1. Various physical situations in stellar atmospheres 
defined by the successive break-down of relations valid in ther- 
modynamic equilibrium. Here the symbol ^ means "not equal 
to" or "not valid" 



TE 


LTE 


NLTE 


Interplanetary 


Iv = By 




Iv^By 


Iv^By 


Boltzmann 


Boltzmann 


T^Boltzmann 


T^Boltzmann 


Saha 


Saha 


T^Saha 


T^Saha 


Maxwell 


Maxwell 


Maxwell 


T^Maxwell 



d. Departures from thermodynamic equilibrium: LTE, NLTE 

Deep in a star one has TE (Table 1). Rising towards the surface, the first concept that 
breaks down is the equality of the intensity and the Planck function, while the Maxwell 
and Boltzmann distributions as well as the Saha equation are still valid. This is called 
local thermodynamic equilibrium (LTE). LTE holds roughly up to the photosphere. 
Rising further into the chromosphere and inner corona, both the Boltzmann distribution 
and the Saha equation are no longer valid. This situation is called Non-LTE (NLTE). 
Here the individual transition rates between the various energy levels must be considered 
in detail. In the low density parts of the corona the temperatures of the different particle 
species become unequal, and eventually in the interplanetary medium even the Maxwell 
distribution breaks down. 



Absorption & Emission Coefficients, Source Function 



Iv 















As 



Iv+AIv 









\ 





FIGURE 6. Left: Absorption coefficient. Right: Emission coefficient 

a. Absorption coefficient 

Consider in Fig. 6, left, a ray of light with intensity /y, penetrating a box with surface area 
AA and thickness /Ss. The absorbing or scattering atoms in the box have a number density 
n [ cm^^ ] and cross section qy [ cm^ ]. The sum of all cross sections is Qv = L^v, thus: 



Iv 



AA 



Y^qy qy n l!A l!^ 



AA 



AA 



-Ky Ai' 



(41) 
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where 

Kv = qv n (42) 

is the absorption coefficient or opacity and has the dimension [cm~^]. The function 
Ky{T,p) is available as a program package, e.g. in the program MULTI discussed below. 
It is often sufficient to consider only a gray (i.e., frequency averaged) Rosseland opacity: 



K = 



/ CO 



CO 

\ 



(43) 



At photospheric and low chromospheric temperatures the gray opacity k can be approx- 
imated by the H~ contribution: 



- = 1.376- lO-V'^^Sj^s 
P 

b. Emission coefficient and source function 



2 —1 

cm g 



(44) 



Let the volume AV (Fig. 6, right) emit photons of energy Afi in direction e„ into the 
solid angle Afl per frequency interval Av and per time interval A?. Then 



77v (r,e„,0 = lim- 



AE 



(45) 

Ay A? An Av 

is the emission coefficient with the dimension [ erg cm~^ s~^ sr~^ Hz~^ ]. Here the limit 
is taken for AV, Af, /S£i, Av — > 0. The source function is defined as 



JCv(r,0 



(46) 



which has the same dimension as the intensity. In TE the amount of energy absorbed, 
AEa, is exactly equal to the amount of energy emitted, AEe 

AEa = By Kv As AA At AO. Av ; AEe = rjv AA As At AO. Av , 

which leads to Kirchhoff's law: 

Sy=By . (47) 

While in NLTE Kirchhoff's law no longer holds, it does hold in LTE because Sy, as we 
will see below, is essentially the population ratio, n„/ n/, and from our definition of LTE, 
this ratio, same as in TE, obeys the Boltzmann distribution. 



Radiative Transfer 

a. Transfer equation 

Consider a gas layer in a stellar atmosphere (Fig. 7). Let e„ and the geometrical height x 
point in the outward vertical direction. Consider a light ray in an arbitrary direction e^. 
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The angle between e„ and e', is Let s be the geometrical distance along this light ray. 
From Eqs. (41), (45), and (46) we have 

d/v = -/v Kv ds + riv ds = - Ky (ly -Sy) ds (48) 

and 

dx= -d^ cos(180°-i?) =d^ cosi^ =ju d^ . (49) 



stellar 
interior 



Stellar 
surface 



ds 




dx 



FIGURE 7. Transfer equation 
This gives the radiative transfer equation: 

dly 



— = - KTv (Iv-Sy) 



We define the optical depth Ty-. 



Ky dx or dTv = — Ky dx . 



Using Eq. (51), the transfer equation can be written in terms of the optical depth: 

^ — = Iv-Sy . 
dZy 

b. Formal solution of the transfer equation 

Multiplying Eq. (52) with e^'^"^^ one finds 



d (lye-^'/^ 
dTv/M 



-Sye 



'Ty/fl 



which can be integrated between the limits a and b: 

b 



lye-'^l^ = - [Sy e-^v/^ d</ju . 

a J 



(50) 



(51) 



(52) 



(53) 



(54) 
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Boundary conditions: 

i) outer boundary: No incoming radiation from outside the star 

/v(0,Ai)=0 , ju<0 , (55) 

ii) inner boundary: At the center of the star the outgoing intensity is finite 

/v (Tvc=o,ju) = finite , jU > . (56) 
Consider the case jU < 0, take b= Ty and a = 0, and multiply Eq. (54) with e^*'/'^: 
ly (Tv,M) e-^'^^ e^'/^" -ly (0,m) e^'/^ = 

ly (Tv, jU) = - / Sy «) ^-«-^v)/M d</;U , (57) 


This is valid for incoming radiation, where ji <0 and 13- > n/2. 

Now consider the other case jU > 0, take b = Zyoo very large and a = Zy, and multiply 
Eq. (54) with-e^v/^: 

-ly (Tvoo,ju)£1^ e'^/^' +Iy (Tv,ju) e-^'/^" e^^/^ = 
=0 

oo 

7v(Tv,ju) = / Sy{z'y) d</;u , (58) 

Tv 

This is valid for outgoing radiation, where jU > and 'd- <7t/2. 

Radiative Equilibrium, Eddington Approximation 

a. Radiative equilibrium 

1 

Operate with lit J djU on Eq. (50) and use Eqs. (32), (35): 
-1 



djc 
where 



" = -An Ky {Jy-Sy) , (59) 



+1 +1 



2n J Sydii = In Sy J dii= An Sy , 



(60) 



-1 -1 
if is assumed independent of jU. 
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If in a plane-parallel time-independent atmosphere the energy transport is by radiation 
only, then F = dv = const, i.e. there is a constant energy flux through all layers. 

This condition is called radiative equilibrium: 



dip f dF r 
— = / —^dv = ~4n / Ky {Jy-Sv)dv = 
dx J dx J 



(61) 





which for gray opacity Ky = k reduces to 7 = 5 as 7c can be taken out of the integration. 



Iv 






Interior Surface 
FIGURE 8. Isotropy of the radiation field as function of depth in the star 

b. Eddington approximation 

Consider the behavior of the intensity with depth (Fig. 8). If in deeper layers the intensity 
is approximately isotropic we get for the Eddington K-integral: 



1 1 1 

isTv = - y Iv^^dii ^ -Iv J M^d^ = ^Iv 



(62) 



-1 -1 
Similarly for the mean intensity at the same depth we obtain 

+1 ^ +1 

A.. T / _ 



1 1 



(63) 



-1 



-1 



We thus find the Eddington approximation: 



1 



Ky--Jv , (64) 
which gives relatively good results even in situations where ly is fairly anisotropic. 



Gray Radiative Equilibrium Atmosphere 

Assume LTE, the Eddington approximation, radiative equilibrium, and gray opacity; 
and integrate all equations over frequency v. Then 



S = B , K=^J , 7 = 5 
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from which we find with Eq. (37): 

j = S = B=-T^ . (66) 

, +1 

Operate with ^ / ju dju on the transfer eq. (52) and use Eqs. (35), (64): 
-1 

+1 +1 
dK \6J \ f If F 

-1 -1 

=F/4k =0 

Multiply with 3 and integrate over t: 

/= - r^= — T + 7(0) . (68) 

At the surface the ingoing intensity is zero, so we need to integrate only over jU > in 
Eqs. (32) and (35): 

1 1 
J{Qi) = \j B{Q,)dn = ^ , F{0) = 27tjB{0)^d^^nB{0) , (69) 



thus 

7(0) = ^ • (70) 

Let us define: 

F = F{0)^aT^,, , (71) 

where T^^ is the effective temperature. Then we find from Eqs. (68), (70), and (71) the 
so called gray radiative equilibrium T{t) relation: 

n4 ^ rr4 lr= 2 



7- =47;ff + . (72) 

The relation between the optical and geometrical depths is obtained from Eq. (51): 

df=-Kdx . (73) 

For a plane, static (where the flow velocity v = 0) atmosphere we find from Eq. (9) the 
equation of hydrostatic equilibrium: 

dp = -p g dx . (74) 

Using for example the simple opacity law (44), it is seen that Eqs. (72), (73), and (74) can 
be integrated as functions of x. This allows to construct a static, plane, gray, radiative 
equilibrium atmosphere if the two quantities T^ff and g as well as suitable boundary 
conditions are given. Together with Eq. (12), p = p%T / pL, we have four equations for 
the four unknowns T,p,p, and T. Such radiative equilibrium atmospheres are the starting 
atmosphere models for chromospheric wave calculations. 
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Net Radiative Heating Rate, Radiative Heating Function 

The net radiative heating rate, 4>r [ erg cm~^ ], (see Eq. 1 1) is the negative diver- 
gence of the total radiative energy flux F = JgFy dv, thus in the ID case we have 

dF 

^R = --r • (75) 
ax 

Integrating Eq. (59) over v we find from Eq. (75) the radiative heating rate 

poo 

^^ = 471 JCv {Jv - Sv) dv , (76) 
^0 



NLTE THERMODYNAMICS AND RADIATION 



As noted above (cf. Table 1), NLTE conditions prevail in the chromosphere and corona, 
so that in these outer stellar regions the individual transitions giving rise to lines and 
continua have to be considered. Chromospheres and coronae are regions where there is 
mechanical heating and where departures from LTE are important. Starting with a slight 
departure in the upper photosphere, NLTE becomes extreme in the high chromosphere 
and corona. It is therefore necessary to review the transition rates and outline the 
methods to treat the thermodynamics and radiation under NLTE conditions. General 
literature for this section includes [25, 26, 27, 29, 30, 31, 32, 33]. 



Transition Rates for Lines and Continua 



Under NLTE conditions, there exists a well-defined kinetic temperature T that is 
determined by the Maxwell velocity distribution, but the Boltzmann distribution and the 
Saha equation (Eqs. 38, 39) are no longer valid. Conservation of particles requires that 
the population rim of level m obeys the time-dependent statistical rate equation 

l~ ^ ■ ~ ^ i^jfn Pm j i C^^) 

where Pab denotes the transition rates (= the number of transitions per sec) from level 
a to level b. In cases where populations adjust faster than the time scale over which 
varies in hydrodynamic changes, the statistical equilibrium equation holds, 

fT-jPjm — ^j{Rjm'\'Cjrn} — nm Pmj — {Rmj '\' Cmj) i (^8) 

where R denotes radiative and C coUisional transition rates [cm~^s~^]. Here i?| is 
the absorption rate, P}^^ the induced emission rate, R^^ the spontaneous emission rate, 
C| the coUisional excitation or ionization rate, and Cj^ the coUisional deexcitation or 
recombination rate. 
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a. Lines 



The radiative transition rates between two bound energy levels, a lower level / and an 
upper level u (Fig. 5), are 



ind 



=niRiu = niBiuJiu , =nuBuiJiu , R{=nuAui 



(79) 



where A„/, and B^i are the Einstein coefficients (tabulated e.g. by Allen [29, 30]). 
Jlu is the mean intensity Jy, averaged over the line. If is the line profile, then 



J 



lu 



(pyjydv with / 9vdv = 1 
Av ./Av 



(80) 



Here the frequency integrals extend over the width Av of the line. In the above equations 
we have assumed that the emission and absorption profiles of the line are identical, that 
is, we assume complete redistribution, CRD. The line profile is usually given by the 
Voigt profile (see Fig. 9) 

ipy = —}—H{a,y) , (81) 
0rAvD 

where the damping parameter a and the normalized frequency separation v are given by 

r V- Vo 



4;rAvD 



Avd 



(82) 



Here Vq is the line center frequency, Avd the Doppler width and T the damping constant. 
In thermal equilibrium (TE), detailed balancing is valid, i.e. the radiation and collision 




3Av^ 

FIGURE 9. The Voigt spectral line profile consists of folded Gauss and Lorentz profiles, 
processes individually balance each other, 

R.^=R'f+Rj , niBiuBy=n,,Aui+n,,BuiBy . 

Dividing by nuBui, solving for By, and using Eqs. (36) and (38), one obtains 

Aui 1 A,,, 1 2hv^ 1 



B, 



g, l^^lu \ B, 



C ekT — 1 



(83) 



(84) 



and because T is arbitrary 



htl 



-B 



id 
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giBiu = guBui . (86) 

From Eq. (79) one sees that A„; has the dimension [ s~^ ], while B^i and 5/„ have the 
dimension [ cm^ sr Hz erg~ ^ ] . Eqs. (85) and (86) are relations between probabilities that 
involve only atomic parameters and do not depend on atomic level populations, they 
thus are also valid in NLTE. Let us introduce an absorption cross section [ cm^ ] 

Ttc^ hvi 
a/„(v) = fi(p^=Biu-^(py . (87) 

Here e is the elementary charge, ntg the electron mass, and // the oscillator strength. Let 
us write the Boltzmann distribution (38) in the form 



4 = ^.-^ , (88) 
SI 

where gi = 2i^ and populations marked by a * indicate quantities in TE or LTE. One 
finds 

2hv^ tit hv n1 hv 

Aui=Bi^^^e-w , Bui=Biu-^e-w . (89) 

With this the radiative transition rates can be written 

[ An 

R] = ni aiu{v)—Jvdv = niRiu, (90) 

JAv "V 



^ nt JAv hv 



n„ JAv 



(92) 



Defining 



Gui = - , Rli= aUv)—[^+Jv]e-^dv , (93) 



lAv hv 

the total radiative deexitation rate is given by 
i?^=7?;P+7?f ^nA, = n„_^^a;„(v)G„,^(^^+7v)dv=n,|i?;, . (94) 

For the collisional excitation and deexcitation rates [ cm"^ s""^ ] one has, respectively, 

C| = ni Ciu = ni n^Q.iu{T) , C| = n„ C„; = n„ n^Q.ui{T) . (95) 

The fls are called collision cross sections. In TE, because of detailed balancing, one 
finds 

Cul = -^Ciu = -^n,ai,{T) = ^e^n^a^T) . (96) 
nu* riu* gu 

This relation between the collision cross sections is also valid in NLTE because it 
involves only atomic parameters and the Maxwell velocity distribution. 
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b. Continua 

The radiative transition rates [cm^-^ s^^ ] between a bound level / and the continuum k 
(Fig. 5) can be derived using similar arguments and detailed balancing in each frequency 
interval: 

a/(v)— iydv , 



'V/ 

4k 



(97) 
(98) 



In TE and LTE the Saha equation (39) is written 

^ \lKmJiT J 2Mjt ' 

where Ei = hVi is the energy difference between level / and the continuum, in Eq. 
(100) is the partition function 

Uk=2^8ie , (101) 

i=l 

where the summation is carried out over the imax bound levels of the next higher 
ionization stage. It can be shown that Eqs. (98) and (99) are also valid in NLTE. Defining 



(102) 



< _ftv _+ r , . 4k flhv'^ 

the total radiative recombination rate is given by 

/"^ 4-k / 2/zv \ /I* 

ai{v)Gki-^i^-^+Jyjdv = nk-^Rli . (103) 

For the collisional ionization and recombination rates [ cm~^ s~^ ] one finds similarly 

C| = niQk = nin^Q-ik (T) , Q = nkCu = nkn^Q-ki (T) = nkn^Q.ki {T, n^) . (104) 
In TE, because of detailed balancing, one has 



n* 



q = UkCki = Uk-^Cik = tik-^n^Q-ik (T) , (105) 



with n*/«| given by Eq. (100). This relation is also valid in NLTE. Partition functions, 
absorption and collision cross sections (a's and Cl's) can be found in Allen [29, 30] and 
references cited there. 

Chromospheric Dynamics and Line Formation May 25, 2008 22 



Line and Continuum Source Functions 



We now consider the transfer of radiation through a stellar gas. After Eq. (46) the 
source function is defined as the ratio of the emission and absorption coefficients. 

a. Lines 

From Eqs. (45) and (79), noting that the dimension of the transition rates is [ cm~^ s~^ ], 
one finds for the line emission coefficient [ erg cm~^ s~^ sr~^ Hz~^ ] 



dRf hv hV 2hV^ 

=^TZ= ^uAuij-(Pv = -^aUv)n,G,i. (106) 



Similarly from Eqs. (41) and (79), noting that 



^=-z...=- , (107, 

As dv 47t 



the line opacity [cm ^ ] is given by 



= {niBiu-nuBui) ^<Pv = {l - ^<Pv = a/«(v)(«/ -n„G„/) . (108) 

The absorption coefficient is defined for a total intensity change (which results from 
absorption minus induced emission). The line source function (see Eq. (46)) is 



5line= ^ nA^i ^ 2hV^ 1 ^ 2hV^ 1 
-^ijne niBiu-n^B^i c2 ^ - 1 c^ ke^-i 



where we have used the Boltzmann distribution. The departure from LTE coefficient bi 
is defined by 

b,^T^ . (110) 

Note that always = 1- It is seen that in LTE 

5}r = 5v, with bi^K^l . (Ill) 
In deriving Eqs. (106) and (108) we have assumed complete redistribution (CRD), 
b. Continua 

From Eqs. (45) and (99) we have similarly a continuum emission coefficient [ erg cm~^ 
s-i sr-i Hz-i ] 



cont = ^ „^ J^oj^ ^ -^e-w = ai{v)nkGki . (112) 
dv An nl c^ c^ 
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Also, using Eqs. (97), (98), and (107), we obtain the continuum opacity [cm ^ ] 

JC^°"t = ai (v) (ni-uk'^e-^^ = a/(v)(n; -n^G^;) . (113) 
The continuum source function is then given by 



„cont 2hV^ 1 2/7 1 

ocont _ '/v _ ^ _ ^ /I i^x 



Note that in LTE 

5^°°^ = 5v, with &/ = l . (115) 

Eqs. (Ill) and (1 15) show the validity of Kirchhoff 's law (47) in LTE, as stated already 
above. Taking the source function equal to the Planck function is used as primary 
definition of LTE in some texts. For us here, the equality of the source and Planck 
functions (Eqs. (1 1 1) and (115)) is the result of the way by which we defined LTE. 



Computation of the LTE and NLTE Level Populations 

Suppose the temperature T and gas pressure p are given along with the chemical ele- 
ment abundance of the stellar gas. How can the number densities and level populations 
of the different atoms and ions be computed? Let 

'max ^max 

nr,s,i , «*,(=J2"W ' «/=12"*'' ' ^^^^^ 

r=l s=l 

be the number density [ cm^^ ] of particles in energy level r and ionization stage s of 
element i, the number density of particles of element / in ionization stage s, and the total 
number density of particles of element /, respectively. Here r = 1 , • • • , rmax. where rmax is 
the number of bound levels, and s= 1 , • • • , 5max> where Smax is the number of ionization 
stages. 

a. LTE populations 

We first assume LTE. In this case the Boltzmann distributions and Saha equations are 
valid. Summing Eq. (38) over the bound levels r we obtain the Boltzmann distributions 
with the partition functions 

Jl£ = ^Jl±e-^ , u^,^Y.^r,s,ie ^ . (117) 

Similarly summing Eq. (39) over all bound levels / we find the Saha equations 

-^±11 = ?±hLl — ^ — e w — . (118) 

ns,i rie Us,i \ J 
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One now proceeds as follows: 1. estimate rie. 2. Compute the element densities = 
{p/kT — n^Ail E;^;- 3. Compute the risX, for this we have ^max — 1 Eqs. (118) and the 
third of the Eqs. (116). Similarly compute rir^s.i using the rmax Eqs. (117). 4. Evaluate the 
new electron density using the equation of charge conservation 

n^ = Y^Y^^^s,i ■ (119) 

Going back to step 1, we could iteratively improve until a converged solution is ob- 
tained. One actually uses the much faster Newton-Raphson method by writing Eq. (119) 
as /(we) = 0. If the true solution is written = «e + ^"e, where /ig is an estimate, 
we have /(^e) + ^^"e = 0. The derivative can be evaluated analytically. Solving for 
5ne = — /("e) Iw^^^ ^ estimate + bn^, etc. This Newton-Raphson iteration 
converges very fast if the initial estimate is reasonable. 

b. NLTE populations 

In NLTE the situation is different. Here not only the Boltzmann distributions and Saha 
equations are no longer valid, but also the mean intensities 7v are unknown. Of the dif- 
ferent methods to solve the problem we discuss only the complete linearization method. 
A computer code for this method, MULTI (Scharmer and Carlsson [31], Carlsson [34], 
Carlsson [32]) is available [33] and widely used in the astrophysical community. 

Here one also has a given temperature T- and pressure p-distribution. One starts by 
selecting an electron density We-distribution and assuming level populations n,-. Then 
the radiation fields Jy and the statistical rate equations are computed, which leads to 

improved n,. In a Newton-Raphson scheme a system of equations, with a matrix W 
operating on the vector of variations { } being equal to an error vector {£■, }, is inverted 
and the new estimates Ui + dui evaluated. After convergence of this iteration the Re- 
distribution is modified until the given p-distribution is obtained. Note that physical 
vectors are directed quantities in space, mathematical vectors are simply arrays denoted 
by curly brackets { }. 

Assume that for the step of the iteration scheme we have the populations n-"^ and 

the transition processes pj-""^ and seek small corrections drii and dPij. From Eq. (78), 
written for level /, with the total number of levels NL = N+l,'we have 

(.W + 5.f))i:(/^W + 5/^W)-£(nW + ^"))(pjf^^ . (120) 

After expanding and neglecting second order terms we get 

Nl Nl Nl Nl 

5nf^ £4") +nW £ SP.f - t^ Snfpf - f^nf 5pf = E^"^ , (121) 

(n) 

where Ej^ ' represents the zeroth order terms and vanishes for the converged solution. 
Our aim is to write dPij, 5Pji in terms of 5ni, and then to solve Eq. (121) for 5ni. With 
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FIGURE 10. Matrix equation in the complete linearization method as used by Carlsson [34] 



5Pij — 5Rij + 5Cij, where SQj — and SR^'^ = because and T are given, we have 



Iff Alt 

5Pij = 5Rij = 2j J ^"O^'/^^vAidvdAi , i> j 

-lAv 

-1 



\l I ^^O^^v/idvdjU , i<j 



(122) 



lAv 



Uij is given by Eq. (87) for lines, and a,^ = a, (v) from Eq. (97) for continua. The Gij 
are given by Eqs. (93) and (102), and the v-integration interval Av is either over the line 
width or from v,- to infinity, depending on whether a line or a continuum transition is 
considered. 

The intensities are obtained from a linearization of the radiative transfer equa- 
tions, where it is important that the in- and outgoing intensities at depth d originate 
from other points of the atmosphere than depth d, depending on the considered fre- 
quency. This is done using a linear system in matrix form shown in Fig. 10. With vectors 
dn = {5nid}, E = {Ej^}, where z = 1, ■ ■ • , A^/^ is the energy level index and J = 1, ■ ■ ■ ,Nd 

the depth index, Eq. (121) can be written with a grand matrix W, 



W dn = E . 



(123) 



As an example Fig. 1 1 shows the non- vanishing elements of the grand matrix V7 for a 5 
level -I- continuum Ca 11 calculation. The non-zero matrix elements are clustered along 
the main diagonal, showing that the most important radiative interactions are of inter- 
mediate range. The off-diagonal matrix elements represent the non-local contributions. 
This band structure permits an efficient matrix inversion. 
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FIGURE 11. Grand matrix W for a 5 level + continuum Ca II calculation of the solar VAL81-C model 
[1], after Carlsson [34]. 



Chromospheric Radiation Loss 

From the procedure discussed above the net heating rates [ergcm^^ s^^ ] for lines 
and continua can be determined as 4>/„, 4>/ = hv — : 



4>/„ = n/4;r / a/„(v)7vdv -n„^4;r / a/„(v) l^^ + Jv] e ^dv 



Av 



Av 



4>/ = n/4;r/ a/(v)7vdv -n^^4;r / ai{v) (^^^^ + Jv] e ^dv 



Vl 



V/ 



(124) 
(125) 



Coronal Radiation Loss 

In the tenuous coronal layers the temperature is high and the density is very low. 
Consider the energy levels of a typical multiply ionized coronal ion with its large 
ionization Ei and excitation energies (Fig. 5). Because the photospheric radiation 
field Jy does not provide photons with enough energy to excite a coronal ion, ^ 
and ~ 0. In addition, as is very small in the corona, C| ~ 0. What remains from 
Eq. (78) is 

C^=RT , (126) 



which is called the thin plasma approximation. From Eqs. (99) and (104): 



nin^ai{T)=nk^ tti v ^dv = nknJ{T) 

nj. Jvi nv 



(127) 



where we used the Saha eq. (100), 



nl _ 1 f2KmM\^'^'^Uk 



(128) 
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FIGURE 12. Ionization ratios for Fe in the thin plasma approximation valid for the corona 



and collected the temperature-dependent parts in some function f{T). This gives 



ni f{T) 



g{T) . (129) 



The function g{T) can be tabulated (e.g. Jain and Narain [35]). Note that as lower bound 
level we took 1 = 1 because only the ground level is significantly populated. Eq. (129) 
states that very differently from the Saha equation (128) in LTE, which depends also on 
He, the coronal ionization ratio depends only on T, due to the thin plasma approximation. 
To illustrate this dependence, Fig. 1 2 shows the ionization ratios ,/ Y,s ^s,i for Fe in the 
coronal approximation. 

The coronal radiation loss due to lines is 

-^YL = ^7n]-AKKl=Y,hvn,,R'^^=Y,hvnionn^Q.{T) , (130) 

where the sum is taken over all lines. The total cooling rate includes also other processes 
and is given by 

-4>R = nH«e/^ad(r) , (131) 

where PKa.d{T) is a function shown in Fig. 13. Similar functions were computed by 
numerous authors (e.g., [36, 37, 38]). The descending part is often approximated by a 
simple power law. A dependence o^T^^/^, like 

fliad- 5- 10-20^-1/2 (132) 

(after [37]) is particularly useful as it allows to cast the coronal energy balance in dimen- 
sionless form, thus making it possible to scale corona models from one star to another 
(Hammer [39]). Please note that Eq. (131) defines only a radiative cooling function, 
which will ultimately cool the gas down to temperatures T ^ in the absence of heat- 
ing. In reality, if T becomes small enough, radiative heating must also be considered, 
similarly as in 4>r = 4;rfc (7 — 5) . 
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FIGURE 13. Radiation function PRad as function of temperature and various numerical fits of this 
function, after McWhirter et al. [37]. 



SIMULATIONS VS. OBSERVATIONS 



General references on the concepts discussed in this chapter include [40, 41, 42, 43]. 



Numerical Simulations 

In order to investigate the dynamics of the chromosphere theoretically, it is necessary 
to solve the conservation equations for mass (18), momentum (20), and energy (25) as 
functions of space and time for given initial and boundary values. We also need the equa- 
tion of state and other thermodynamic relations necessary to close the system, Eqs. (12) - 
(14). Simultaneously we must calculate the ionization rates, the electron density, and the 
level population densities and radiative energy source and sink terms, Eqs. (124) and 
(125), for all important spectral lines and continua. As initial atmosphere one can e.g. 
use a gray radiative equilibrium atmosphere, as derived above. And as boundary con- 
ditions one typically specifies the velocity at the lower end of the computational region 
and allows for waves to leave the upper end with as little reflection as possible. 

Ideally one would like to solve this system of partial differential equations in three 
dimensions (3D), in order to be able to handle the chromospheric structure. Such 3D 
simulations exist (e.g., Wedemeyer et al. [44]), in some cases even including the effects 
of magnetic fields (S chaff enberger et al. [45], Steiner et al. [46], Hansteen et al. [47]) 
- but unfortunately the currently available computer power does not yet permit the 3D 
treatment of the full problem with all the physics that is important in the chromosphere. 
Radiative transport must either be treated in gray LTE, or if NLTE effects are considered, 
they must be highly simplified; moreover the spatial resolution of fine structures such 
as shock waves is poor. Nevertheless, these 3D simulations provide impressive results 
about the structuring and dynamics of the chromosphere (see also Steiner, this volume). 

Another approach restricts itself to one spatial dimension (ID), but solves the full 
NLTE radiation-hydrodynamics problem and calculates the variation of line profiles. 
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which can then be compared with observations. The two most sophisticated numerical 
codes to deal with this problem are the one developed by M. Carlsson and R. F. Stein 
(briefly described in [43, 48, 49, 50]) and the one developed by P. Ulmschneider and 
collaborators already since the late 1970s (e.g., [51]). We will summarize the basic 
working of the latter code, the most recent version of which is explained in detail in 
Rammacher and Ulmschneider [41]. 

In order to solve a system of partial differential equations, one basically replaces 
differentials by finite differences that connect quantities of the known solution at the 
previous time level with those at the new time level, and then one solves for the latter in 
order to advance the solution in time. 

A special feature of the Ulmschneider-Rammacher approach is that before this is 
done the differential equations are first transformed into their characteristic form; i.e. 
one takes explicitly into account that matter and hydrodynamic signals travel with 
speeds v and v ± c, respectively, thus defining the so-called characteristics in space- 
time (cf. [17, 51]). The method of characteristics has a number of advantages: It is 
numerically efficient [52], and it makes it easy to recognize the formation of shocks 
(namely, when two neighboring characteristics of the same kind intersect) and to handle 
the jump conditions at the shock exactly - i.e., to take care of the continuity of mass, 
momentum, and energy fluxes across the shock, while other variables are allowed to 
change discontinuously. Most other methods cannot treat these discontinuities and use 
artificial or numerical viscosity to spread out shocks over a certain height range, and then 
represent each shock by a number of narrowly spaced grid points. For these reasons, 
characteristics methods with detailed shock handling are very fast. They are, however, 
best suited for ID calculations, since the necessary bookkeeping of characteristics gets 
prohibitively complicated in 3D. 

In this particular code, the solution is advanced in time in an iterative process. Suppose 
that all variables are known at some time level t. At the first time step these are the 
initial conditions to be specified. The radiative heating rates are first assumed to remain 
constant, and a first guess is used for the values of the hydrodynamic variables at a later 
time t + At. This allows to calculate the characteristics and the change of the variables 
along these characteristics, leading to improved values of the hydrodynamic variables 
in the next iteration. This hydrodynamic iteration converges after a few steps. After 
the hydrodynamic variables at the new time level are known, the population levels and 
radiative intensities can be calculated, providing the radiative heating rates. These are 
used in the next hydrodynamic iteration to further improve the hydrodynamic variables, 
which in turn lead to improved population levels and radiative heating rates, and so on. If 
the radiative terms or any other quantities are found to change too rapidly or too slowly 
during this process, the time step A? is decreased or increased accordingly. 

For reasons of computing time efficiency, the radiative part is usually simplified 
during the simulation, and full line profiles for diagnostic purposes are calculated with 
codes like MULTI only at specific time steps. 

A movie clip from an example calculation (Rammacher [53]) shows how a spectrum 
of waves moves through the atmosphere, steepening into shocks that continue to grow, 
whereby occasionally larger ones catch up and merge with smaller ones. The associated 
emission is rather complex: Even though a major contribution comes usually from 
behind large shocks, because of the long-range interaction of chromospheric radiation 
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other parts of the atmosphere can also make significant contributions, depending on the 
changing availability of emitters/absorbers and photons in various parts of the highly 
dynamic atmosphere. In the movie clip [53] this is illustrated for two wavelengths in the 
blue wing of Ca ii K. 



Comparison with Observations 

Such ID simulations have been very successful (Rammacher and Ulmschneider [54], 
Carlsson and Stein [49, 55]) in explaining the Ca H and K Bright Grain phenomenon, 
a characteristic variation of the line profiles of Ca ll H and K that arises when strong 
acoustic shocks traverse the mid-chromosphere. 

Simulations in which the waves were injected according to photospheric velocity 
measurements led to the surprising result that the average temperature beyond a height 
of 500 km (the canonical location of the temperature minimum) did not increase, but 
continued to drop outward (Carlsson and Stein [43, 48]). The emission behind strong 
shocks was found to be so large that no general temperature rise was needed to generate 
the chromospheric emission in lines such as Ca ll H and K. There has been some debate 
if this low average temperature is real or caused by the neglect of high frequency waves, 
which are on principle not observable due to their short wavelengths (Kalkofen et al. 
[56], Carlsson [40]). The problem could also be related to the 3D character of wave 
propagation in the real solar chromosphere, which according to Ulmschneider et al. 
[57] reduces the generation of very strong shocks by the merging of weaker ones, as 
commonly found in ID models, but rather leads to a more continuous heating by a 
larger number of weak shocks. Moreover, the meaning of an "average" temperature 
becomes questionable in the presence of large temperature fluctuations (Carlsson and 
Stein [48], Rammacher and Cuntz [58]). 

Acoustic shock waves have long [59] been thought to be the main heating agent of 
outer stellar atmospheres, or at least of the nonmagnetic parts of stellar chromospheres 
(for a review see e.g. Ulmschneider and Musielak [8]). This view has been challenged 
recently, when Fossum and Carlsson [60, 61, 62] used ID simulations to interpret the 
fluctuations observed by the TRACE satellite in UV continua formed in the upper pho- 
tosphere. They concluded that the small observed fluctuations permit only an acoustic 
energy flux at least an order of magnitude too small to balance the chromospheric en- 
ergy losses. Wedemeyer-Bohm et al. [63] demonstrated, however, that this is mostly due 
to the fact that the spatial resolution of TRACE is insufficient to resolve the fine-scaled 
lateral structuring found in 3D simulations. And Cuntz et al. [64] argued that such a 
low acoustic energy flux would be inconsistent with theoretical calculations of sound 
generation (Musielak et al. [65]) and with the excellent agreement between simulations 
and the measured Ca emission of the most inactive stars. Therefore, it appears that the 
acoustic heating theory is still valid and cannot be considered dead at this time [66]. 

Even though ID simulations describe reasonably well the dynamics of mostly ver- 
tically propagating strong shocks (Ca Bright Grains, as discussed above), they fail to 
match the overall dynamics of the solar chromosphere. If the chromosphere were dom- 
inated by upward propagating plane-parallel waves, the resulting perturbations in the 
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FIGURE 14. Ca II H correlation matrices for VTT observations (top panel) and for a ID model simula- 
tion using a power spectrum with total input energy flux of 5.0 x 10^ ergcm^^ s^^ (bottom panel); from 
[19]. 

Ca II H and K lines would always start in the line wings (which are formed in the pho- 
tosphere) and then move inward towards the line center (which is formed in the upper 
chromosphere). Observations [20] show that this may happen indeed, however one often 
sees large portions of the blue or red wing increase simultaneously. To quantify this be- 
havior, W. Rammacher [19] has introduced a diagram that may serve as a "fingerprint" 
of the dynamics of the chromosphere (Fig. 14). It shows the correlation r of the intensity 
at any wavelength with the intensity at all other wavelengths in the line profile. If two 
parts of the atmosphere, where two different wavelengths are formed, always vary si- 
multaneously, the correlation r is 1. Obviously this must be true for the diagonal. On the 
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other hand, if the dynamics in the two parts of the atmosphere are completely unrelated, 
the correlation is r = 0; whereas r — —I means total anticorrelation, i.e., one part of 
the atmosphere always brightens when the other darkens. Observations (upper panel in 
Fig. 14) show a much higher degree of intensity correlation between different parts of 
the atmosphere than ID numerical simulations (lower panel). The difference to the ob- 
served fingerprint turned out to be large for all numerical simulations, irrespective of the 
amount of wave energy and the spectrum of wave frequencies used. (According to [67] 
the chromosphere depends less sensitively on the wave spectrum than on the energy 
flux.) Therefore, vertically propagating shock waves cannot explain the observations. 
Rammacher et al. [19] suggest that oblique shock fronts could explain the high observed 
correlations, because they would lead to a simultaneous brightening of deeper and higher 
parts of the atmosphere. Such oblique shock fronts are commonly found in 3D simula- 
tions (Wedemeyer et al. [44], Wedemeyer-Bohm et al. [68]); also Ulmschneider et al. 
[57] argued theoretically that 3D propagation of shock waves must be important. 

To summarize, the highly dynamic solar chromosphere (described in the first chapter 
of this paper) calls for the simultaneous solution of the radiation hydrodynamic equa- 
tions under NLTE conditions (as discussed in the next three chapters). Modern comput- 
ers allow their full solution in ID, while for 3D calculations, due to the present lack of 
sufficient computational power, simplifications have to be made. Since 3D effects turned 
out to be important in the real solar chromosphere (as discussed in this chapter), we will 
have to take advantage of the best features of both types of calculations over the next 
few years and combine them with the ever improving observational capabilities in order 
to understand the dynamic solar chromosphere. 
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